#    pythonequations is a collection of equations expressed as Python classes
#    Copyright (C) 2008 James R. Phillips
#    2548 Vera Cruz Drive
#    Birmingham, AL 35235 USA
#    email: zunzun@zunzun.com
#
#    License: BSD-style (see LICENSE.txt in main source directory)
#    Version info: $Id: __init__.py 267 2010-09-25 13:25:43Z zunzun.com $

import pythonequations, pythonequations.EquationBaseClasses, pythonequations.ExtraCodeForEquationBaseClasses
import numpy, scipy.interpolate
numpy.seterr(all = 'raise') # numpy raises warnings, convert to exceptions to trap them


class Spline3D(pythonequations.EquationBaseClasses.Equation3D):
    RequiresAutoGeneratedGrowthAndDecayForms = False
    RequiresAutoGeneratedOffsetForm = False
    RequiresAutoGeneratedReciprocalForm = False
    RequiresAutoGeneratedInverseForms = False
    _name ="Spline"
    _HTML = "z = B-Spline Interpolation Surface"
    splineFlag = 1
    function_cpp_code = ';' # unused

    def CreateCacheGenerationList(self):
        self.CacheGenerationList = []
        self.CacheGenerationList.append([pythonequations.ExtraCodeForEquationBaseClasses.CG_X(NameOrValueFlag=1), []])
        self.CacheGenerationList.append([pythonequations.ExtraCodeForEquationBaseClasses.CG_Y(NameOrValueFlag=1), []])


    def EvaluateCachedData(self, coeff, _id):
        return self.scipySpline.ev(_id[0], _id[1])


    def FitToCacheData(self):
        self.FindOrCreateCache()
        self.scipySpline = scipy.interpolate.fitpack2.SmoothBivariateSpline(self.cache[0], self.cache[1], self.DependentDataArray, s=self.smoothingFactor, kx=self.xOrder, ky=self.yOrder)
        self.coefficientArray = self.scipySpline.get_coeffs() # scipySpline now has spline coeffs internally, also make them available as "coefficientArray"


    def CalculateFittingTarget(self, in_coeffArray):
        raise NotImplementedError, 'Not implemented for spline interpolations'


    def CodePYTHON(self):
        s  = '''# To the best of my knowledge this code is correct.
# If you find any errors or problems please contact
# me at zunzun@zunzun.com.
#      James
#
# the code below was partially based on the fortran code at:
# http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbisp.f
# http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbspl.f

'''
        s += "def " + self.__class__.__name__ + "_model(x_in, y_in):\n"
        s += '    tx = ['
        for i in range(len(self.scipySpline.get_knots()[0])):
            s += "%-.16E" % (self.scipySpline.get_knots()[0][i])
            if i < (len(self.scipySpline.get_knots()[0]) - 1):
                s += ", "
        s += ']\n'
        s += '    ty = ['
        for i in range(len(self.scipySpline.get_knots()[1])):
            s += "%-.16E" % (self.scipySpline.get_knots()[1][i])
            if i < (len(self.scipySpline.get_knots()[1]) - 1):
                s += ", "
        s += ']\n'
        s += '    coeff = ['
        for i in range(len(self.scipySpline.get_coeffs())):
            s += "%-.16E" % (self.scipySpline.get_coeffs()[i])
            if i < (len(self.scipySpline.get_coeffs()) - 1):
                s += ", "
        s += ']\n'
        s += '    nx = ' + str(len(self.scipySpline.get_knots()[0])) + '\n'
        s += '    ny = ' + str(len(self.scipySpline.get_knots()[1])) + '\n'
        s += '    kx = ' + str(self.xOrder) + '\n'
        s += '    ky = ' + str(self.yOrder) + '\n'
        s += '''
    h = [0.0] * 25
    hh = [0.0] * 25
    w_x = [0.0] * 25
    w_y = [0.0] * 25

    kx1 = kx+1
    nkx1 = nx-kx1
    l = kx1
    l1 = l+1

    while x_in >= tx[l1-1] and l != nkx1:
        l = l1
        l1 = l+1
        
    h[0] = 1.0
    for j in range(1, kx+1):
        for i in range(j):
            hh[i] = h[i]
        h[0] = 0.0
        for i in range(j):
            li = l+i
            lj = li-j
            if tx[li] != tx[lj]:
                f = hh[i] / (tx[li] - tx[lj])
                h[i] = h[i] + f * (tx[li] - x_in)
                h[i+1] = f * (x_in - tx[lj])
            else:
                h[i+1-1] = 0.0
                
    lx = l-kx1
    for j in range(kx1):
        w_x[j] = h[j]

    ky1 = ky+1
    nky1 = ny-ky1
    l = ky1
    l1 = l+1

    while y_in >= ty[l1-1] and l != nky1:
        l = l1
        l1 = l+1
        
    h[0] = 1.0
    for j in range(1, ky+1):
        for i in range(j):
            hh[i] = h[i]
        h[0] = 0.0
        for i in range(j):
            li = l+i
            lj = li-j
            if ty[li] != ty[lj]:
                f = hh[i] / (ty[li] - ty[lj])
                h[i] = h[i] + f * (ty[li] - y_in)
                h[i+1] = f * (y_in - ty[lj])
            else:
                h[i+1-1] = 0.0

    ly = l-ky1
    for j in range(ky1):
        w_y[j] = h[j]

    l = lx*nky1
    for i1 in range(kx1):
        h[i1] = w_x[i1]
        
    l1 = l+ly
    temp = 0.0
    for i1 in range(kx1):
        l2 = l1
        for j1 in range(ky1):
            l2 = l2+1
            temp = temp + coeff[l2-1] * h[i1] * w_y[j1]
        l1 = l1+nky1
        
    return temp
'''
        return s


    def CodeCPP(self):
        s  = '''// To the best of my knowledge this code is correct.
// If you find any errors or problems please contact
// me at zunzun@zunzun.com.
//      James
//
// the code below was partially based on the fortran code at:
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbisp.f
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbspl.f

'''
        s += "double " + self.__class__.__name__ + "_evaluation(double x_in, double y_in)\n"
        s += '{\n'
        s += '    double tx [] = {'
        for i in range(len(self.scipySpline.get_knots()[0])):
            s += "%-.16E" % (self.scipySpline.get_knots()[0][i])
            if i < (len(self.scipySpline.get_knots()[0]) - 1):
                s += ", "
        s += '};\n'
        s += '    double ty [] = {'
        for i in range(len(self.scipySpline.get_knots()[1])):
            s += "%-.16E" % (self.scipySpline.get_knots()[1][i])
            if i < (len(self.scipySpline.get_knots()[1]) - 1):
                s += ", "
        s += '};\n'
        s += '    double coeff [] = {'
        for i in range(len(self.scipySpline.get_coeffs())):
            s += "%-.16E" % (self.scipySpline.get_coeffs()[i])
            if i < (len(self.scipySpline.get_coeffs()) - 1):
                s += ", "
        s += '};\n'
        s += '    int nx = ' + str(len(self.scipySpline.get_knots()[0])) + ';\n'
        s += '    int ny = ' + str(len(self.scipySpline.get_knots()[1])) + ';\n'
        s += '    int kx = ' + str(self.xOrder) + ';\n'
        s += '    int ky = ' + str(self.yOrder) + ';\n'
        s += '''
    double h [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    double hh [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    double w_x [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    double w_y [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

    int i, j, li, lj, lx, ky1, nky1, ly, i1, j1, l2;
    double f, temp;

    int kx1 = kx+1;
    int nkx1 = nx-kx1;
    int l = kx1;
    int l1 = l+1;

    while ((x_in >= tx[l1-1]) && (l != nkx1))
    {
        l = l1;
        l1 = l+1;
    }
    
    h[0] = 1.0;
    for (j = 1; j < kx+1; j++)
    {
        for (i = 0; i < j; i++)
        {
            hh[i] = h[i];
        }
        h[0] = 0.0;
        for (i = 0; i < j; i++)
        {
            li = l+i;
            lj = li-j;
            if (tx[li] != tx[lj])
            {
                f = hh[i] / (tx[li] - tx[lj]);
                h[i] = h[i] + f * (tx[li] - x_in);
                h[i+1] = f * (x_in - tx[lj]);
            }
            else
            {
                h[i+1-1] = 0.0;
            }
        }
    }
    
    lx = l-kx1;
    for (j = 0; j < kx1; j++)
    {
        w_x[j] = h[j];
    }

    ky1 = ky+1;
    nky1 = ny-ky1;
    l = ky1;
    l1 = l+1;

    while ((y_in >= ty[l1-1]) && (l != nky1))
    {
        l = l1;
        l1 = l+1;
    }
    
    h[0] = 1.0;
    for (j = 1; j < ky+1; j++)
    {
        for (i = 0; i < j; i++)
        {
            hh[i] = h[i];
        }
        h[0] = 0.0;
        for (i = 0; i < j; i++)
        {
            li = l+i;
            lj = li-j;
            if (ty[li] != ty[lj])
            {
                f = hh[i] / (ty[li] - ty[lj]);
                h[i] = h[i] + f * (ty[li] - y_in);
                h[i+1] = f * (y_in - ty[lj]);
            }
            else
            {
                h[i+1-1] = 0.0;
            }
        }
    }

    ly = l-ky1;
    for (j = 0; j < ky1; j++)
    {
        w_y[j] = h[j];
    }

    l = lx*nky1;
    for (i1 = 0; i1 < kx1; i1++)
    {
        h[i1] = w_x[i1];
    }
        
    l1 = l+ly;
    temp = 0.0;
    for (i1 = 0; i1 < kx1; i1++)
    {
        l2 = l1;
        for (j1 = 0; j1 < ky1; j1++)
        {
            l2 = l2+1;
            temp = temp + coeff[l2-1] * h[i1] * w_y[j1];
        }
        l1 = l1+nky1;
    }
        
    return temp;
}
'''
        return s


    def CodeJAVA(self):
        s  = '''// To the best of my knowledge this code is correct.
// If you find any errors or problems please contact
// me at zunzun@zunzun.com.
//      James
//
// the code below was partially based on the fortran code at:
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbisp.f
// http://svn.scipy.org/svn/scipy/trunk/scipy/interpolate/fitpack/fpbspl.f

'''
        s += "class " + self.__class__.__name__ + "\n"
        s += '{\n'
        s += "    double " + self.__class__.__name__ + "_evaluation(double x_in, double y_in)\n"
        s += '    {\n'
        s += '        double tx [] = {'
        for i in range(len(self.scipySpline.get_knots()[0])):
            s += "%-.16E" % (self.scipySpline.get_knots()[0][i])
            if i < (len(self.scipySpline.get_knots()[0]) - 1):
                s += ", "
        s += '};\n'
        s += '        double ty [] = {'
        for i in range(len(self.scipySpline.get_knots()[1])):
            s += "%-.16E" % (self.scipySpline.get_knots()[1][i])
            if i < (len(self.scipySpline.get_knots()[1]) - 1):
                s += ", "
        s += '};\n'
        s += '        double coeff [] = {'
        for i in range(len(self.scipySpline.get_coeffs())):
            s += "%-.16E" % (self.scipySpline.get_coeffs()[i])
            if i < (len(self.scipySpline.get_coeffs()) - 1):
                s += ", "
        s += '};\n'
        s += '        int nx = ' + str(len(self.scipySpline.get_knots()[0])) + ';\n'
        s += '        int ny = ' + str(len(self.scipySpline.get_knots()[1])) + ';\n'
        s += '        int kx = ' + str(self.xOrder) + ';\n'
        s += '        int ky = ' + str(self.yOrder) + ';\n'
        s += '''
        double h [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        double hh [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        double w_x [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        double w_y [] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

        int i, j, li, lj, lx, ky1, nky1, ly, i1, j1, l2;
        double f, temp;

        int kx1 = kx+1;
        int nkx1 = nx-kx1;
        int l = kx1;
        int l1 = l+1;

        while ((x_in >= tx[l1-1]) && (l != nkx1))
        {
            l = l1;
            l1 = l+1;
        }
        
        h[0] = 1.0;
        for (j = 1; j < kx+1; j++)
        {
            for (i = 0; i < j; i++)
            {
                hh[i] = h[i];
            }
            h[0] = 0.0;
            for (i = 0; i < j; i++)
            {
                li = l+i;
                lj = li-j;
                if (tx[li] != tx[lj])
                {
                    f = hh[i] / (tx[li] - tx[lj]);
                    h[i] = h[i] + f * (tx[li] - x_in);
                    h[i+1] = f * (x_in - tx[lj]);
                }
                else
                {
                    h[i+1-1] = 0.0;
                }
            }
        }
        
        lx = l-kx1;
        for (j = 0; j < kx1; j++)
        {
            w_x[j] = h[j];
        }

        ky1 = ky+1;
        nky1 = ny-ky1;
        l = ky1;
        l1 = l+1;

        while ((y_in >= ty[l1-1]) && (l != nky1))
        {
            l = l1;
            l1 = l+1;
        }
        
        h[0] = 1.0;
        for (j = 1; j < ky+1; j++)
        {
            for (i = 0; i < j; i++)
            {
                hh[i] = h[i];
            }
            h[0] = 0.0;
            for (i = 0; i < j; i++)
            {
                li = l+i;
                lj = li-j;
                if (ty[li] != ty[lj])
                {
                    f = hh[i] / (ty[li] - ty[lj]);
                    h[i] = h[i] + f * (ty[li] - y_in);
                    h[i+1] = f * (y_in - ty[lj]);
                }
                else
                {
                    h[i+1-1] = 0.0;
                }
            }
        }

        ly = l-ky1;
        for (j = 0; j < ky1; j++)
        {
            w_y[j] = h[j];
        }

        l = lx*nky1;
        for (i1 = 0; i1 < kx1; i1++)
        {
            h[i1] = w_x[i1];
        }
            
        l1 = l+ly;
        temp = 0.0;
        for (i1 = 0; i1 < kx1; i1++)
        {
            l2 = l1;
            for (j1 = 0; j1 < ky1; j1++)
            {
                l2 = l2+1;
                temp = temp + coeff[l2-1] * h[i1] * w_y[j1];
            }
            l1 = l1+nky1;
        }
            
        return temp;
    }
}
'''
        return s


    def CodeCS(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'


    def CodeSCILAB(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'


    def CodeMATLAB(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'

    
    def CodeVBA(self):
        raise NotImplementedError, 'Not implemented for spline interpolations'
